chr po Gene Basal Her2 LumA LumB Normal
1 1 212873846 CENPF 0.4829763 -0.02926616 -0.5437402 0.27822856 -0.07058307
2 1 120072157 PHGDH 0.6345189 -0.18662586 -0.3986822 -1.03013932 0.66043775
3 1 43599336 CDC20 0.3996952 0.00835552 -0.4690440 -0.07041247 -0.04550481
4 1 44992051 KIF2C 0.2035726 -0.16510205 -0.5053947 -0.18289071 -0.39001448
5 1 161575262 NUF2 0.4724002 -0.02381921 -0.7125208 0.58962688 -0.37053337
6 1 200572562 UBE2T 0.3899089 0.28453681 -0.5392594 0.73895213 -0.95238101Workshop MOD: Visualisation de Données Omiques avec R
Présentation du workshop
Ce workshop MOD a été effectué par des étudiants du master Statistiques et Sciences des Données :
- Rim Alhajal : rim.alhajal@etu.umontpellier.fr
- Abdoulaye Diop : abdoulaye.diop03@etu.umontpellier.fr
- Pauline Dusfour-Castan : pauline.dusfour-castan@etu.umontpellier.fr
- Lilou Zulewski : lilou.zulewski01@etu.umontpellier.fr
En génomique, les bases de données à étudier sont généralement très grandes : des visualisations exploratoires permettent de les aborder. Cependant, ces visualisations deviennent complexes dû à la dimension ce qui réduit l’accès aux informations significatives, nuisant aux interprétations. L’objectif de ce workshop est d’améliorer vos compétences en visualisation des données appliquées au cancer du sein.
Données
Les données utilisées dans ce Workshop sont disponibles dans la librairie OmicCircos de R. Elles reposent sur les bases de données fournies par The Cancer Genome Atlas Program (TCGA) contenant des données de patients atteints d’un cancer du sein. Voici la liste des tableaux utilisés :
- TCGA.PAM50_genefu_hg18 : Ce tableau contient un ensemble de 50 gènes identifiés comme une signature d’expression génique associée aux sous-types du cancer du sein (Normal, LumA, LumB, Basal ou Her2) en fonction des schémas d’expression génique.
- TCGA.BC.fus : Ce tableau contient les données de fusion de gènes pour un individu précis.
chr1 po1 gene1 chr2 po2 gene2 name
1 2 63456333 WDPCP 10 37493749 ANKRD30A TCGA-BH-A18M-01A
2 18 14563374 PARD6G 21 14995400 POTED TCGA-BH-A1ET-01A
3 10 37521495 ANKRD30A 3 49282645 CCDC36 TCGA-BH-A1ET-01A
4 10 37521495 ANKRD30A 7 100177212 LRCH4 TCGA-BH-A1EU-01A
5 18 18539803 ROCK1 18 112551 PARD6G TCGA-BH-A1EW-11B
6 12 4618159 C12orf4 18 1514414 PARD6G TCGA-C8-A12Q-01ATCGA.BC.cnv.2k.60 : Les variables de ce tableau sont respectivement le chromosome, sa position, le nom du gène observé ainsi que 60 autres variables représentant 60 individus différents.
TCGA.BC.gene.exp.2k.60 : Sur 60 individus, une métrique liée à l’expression génique et le chromosome, sa position et le nom du gène observé sont consignés.
TCGA.BC.sample60 : Ce tableau référence le type de cancer (Normal, LumA, LumB, Basal ou Her2) dont sont atteints 60 personnes.
TCGA.BC_Her2_cnv_exp : Dans ce tableau, les t-value et p-value de tests statistiques et les métriques FDR et Bonferroni sont référencées.
Recueil des données
Afin de mettre en oeuvre les différentes méthodes de visualisation, il nous faut d’abord recueillir les données et les mettre en forme. Pour cela, suivez les étapes présentées ci-dessous :
- ÉTAPE 1 : Chargement de la librairie
OmicCircos
Pour pouvoir utiliser ÒmicCircos, il faut tout d'abord l'avoir installé à partir de [Bioconductor](https://bioconductor-org.translate.goog/packages/release/bioc/html/OmicCircos.html?_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=tc/) puis l'appeler à l'aide de la fonctionlibrary()`.
Installer Bioconductor
if (!require("BiocManager", tranquillement = TRUE))
install.packages("BiocManager")
BiocManager::install("OmicCircos")library(OmicCircos)- ÉTAPE 2 : Importation des données
Nos données provenant du package OmicCircos (maintenant installé), nous les importons à l’aide de la fonction data().
data("TCGA.PAM50_genefu_hg18")
data("TCGA.BC.fus")
data("TCGA.BC.cnv.2k.60")
data("TCGA.BC.gene.exp.2k.60")
data("TCGA.BC.sample60")- ÉTAPE 3 : Mise en forme des données
# vecteur comprenant les positions de lignes des personnes
# possédant un cancer de type Her2
Her2.i = which(TCGA.BC.sample60[,2] == "Her2")
# vecteur comprenant les identifiants des personnes
# présentant un cancer de type Her2
Her2.n = TCGA.BC.sample60[Her2.i,1]
# indices des colonnes de TCGA.BC.cnv.2k.60 dont
# les identifiants sont présents dans Her2.n
Her2.j = which(colnames(TCGA.BC.cnv.2k.60)%in%Her2.n)
# création d'un data d'expression des différents gènes
# pour chaque individus présentant un cancer Her2
cnv = TCGA.BC.cnv.2k.60[,c(1:3,Her2.j)]
# data cnv auquel on a enlevé les 3 premières colonnes
cnv.m = cnv[,c(4:ncol(cnv))]
# bornage des valeurs >2 et <-2
cnv.m[cnv.m > 2] = 2
cnv.m[cnv.m < -2] = -2
# on ajoute aux 3 premières colonnes de cnv nos données bornées
cnv = cbind(cnv[,1:3], cnv.m)
# extraction des colonnes 1 à 3 de TCGA.BC.gene.exp.2k.60
# ainsi que des colonnes spécifiées par Her2.j
gene.exp = TCGA.BC.gene.exp.2k.60[,c(1:3,Her2.j)]Visualisations classiques
L’objectif de ces premières visualisations est d’assimiler la construction de graphes en reprenant les concepts de base de la librairie ggplot2 de R. Nous vous apprendrons notamment à tracer des graphes de densités et des histogrammes.
L’ensemble des visualisations de cette partie s’appuient sur la base de données TCGA.PAM50_genefu_hg18 et nécessite le chargement de divers packages.
# chargement des packages nécessaires
library(ggplot2)
library(dplyr)
library(tidyr)
library(plotly)
library(viridis)
library(viridisLite)
library(hrbrthemes)Les visualisations permettent d’étudier les différents types de cancer présentés dans le lexique.
Utilisation de ggplot2
Le package ggplot2 fonctionne par couches successives. La première d’entre elles, est un peu le canevas du graphe : elle consiste à indiquer dans quel jeu de données se trouve les données, et quelles sont les variables à représenter. Ensuite, une seconde couche est ajoutée, elle consiste à indiquer le type de graphe à réaliser : scatterplot, boxplot, barplot, etc… Enfin, les couches d’affinage permettent de choisir les couleurs, les échelles des axes, les options de légende et toutes autres options visuelles.
Il ne faut pas oublier le + entre les couches successives.
Comparaison de deux densités
Commençons par étudier certaines densités relatives aux types de cancer du sein deux à deux afin d’en retirer leurs différences.
Voici un exemple très basique de visualisation graphique de la comparaison des densités des types de cancer LumA et LumB qui sont ceux les plus communs :
# création de la base de données
data <- data.frame(value=c(TCGA.PAM50_genefu_hg18$LumA,
TCGA.PAM50_genefu_hg18$LumB),
type=c(rep("LumA", 50),
rep("LumB", 50)))
# création du graphe
comparaison <- data %>%
# création de l'objet
ggplot(aes(x=value, fill=type)) +
# définition du type de graphe
geom_density()
# affichage du graphe
print(comparaison)%>%
Cet opérateur est un pipe, fréquemment représenté par une barre verticale | dans les terminaux : il renvoie la sortie d’une commande vers l’entrée d’une autre.
De nombreuses options de ggplot peuvent améliorer l’apparence générale des graphes. Dans l’exemple ci-dessous, vous pouvez retrouver respectivement la définition manuelle des couleurs, l’option de non-affichage de la légende associée à la coloration en fonction du type, ainsi que la personnalisation des éléments de texte et des légendes. Enfin, l’option alpha de geom_density règle la transparence de remplissage sous les courbes de densité.
# création du graphe
comparaison <- data %>%
# création de l'objet
ggplot(aes(x=value, fill=type)) +
# définition du type de graphe
geom_density(aes(color=type), alpha=0.4) +
# réglage des couleurs de légende
scale_color_manual(values=c("darkblue", "pink")) +
scale_fill_manual(values=c("darkblue", "pink")) +
# choix de l'affichage de la légende
guides(color="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("ARNm expression") +
ylab("Density") +
labs(title="Comparison of gene expression density across different types of cancer",
fill="Cancer type")
# affichage du graphe
print(comparaison)
À l’aide de ce graphe, nous pouvons remarquer que l’ARNmessager des gènes des deux types de cancer LumA et LumB ne s’exprime pas de la même façon. En effet, celui de LumA présente un mode autour de -0.5 tandis que celui de LumB présente un mode légèrement plus haut autour de 0.4.
Il est également possible de représenter les densités sur deux graphes différents en utilisant la commande facet_wrap(~XXX).
Utilisation de
facet_wrap()
Les échelles des différents graphiques peuvent être formatées différemment à l’aide de l’option scales : "fixed" permet d’obtenir des échelles égales pour tous les graphes mais d’autres options telles que "free_x", "free_y" ou "free" sont envisageables.
# création du graphe
comparaison <- data %>%
# création de l'objet
ggplot(aes(x=value, fill=type)) +
# définition du type de graphe
geom_density(aes(color=type), alpha=0.4) +
# séparation du graphe en sous parties
facet_wrap(~type, scales="fixed") +
# réglage des couleurs de légende
scale_color_manual(values=c("darkblue", "pink")) +
scale_fill_manual(values=c("darkblue", "pink")) +
# choix de l'affichage de la légende
guides(color="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size = 8),
plot.title=element_text(size=10, hjust = 0.5, face="bold")) +
# titres
xlab("ARNm expression") +
ylab("Density") +
labs(title="Comparison of gene expression density across different types of cancer",
fill="Cancer type")
# affichage du graphe
print(comparaison)
À vous de jouer
En vous appuyant sur l’exemple précédent, remplissez les trous du code ci-dessous afin de représenter graphiquement la comparaison des densités d’expression d’ARNmessager pour les gènes des types de cancer Her2 et Normal.
Les données relatives aux cancer de type Normal sont contenues dans la variable Normal de la base de données.
Par la suite, les couleurs utilisées pour la représentation de Her2 et Normal seront respectivement darkred et darksalmon.
# création de la base de données
data <- data.frame(value=c(TCGA.PAM50_genefu_hg18$XXX,
TCGA.PAM50_genefu_hg18$XXX),
type=c(rep("XXX", 50),
rep("XXX", 50)))
# création du graphe
comparaison <- data %>%
# création de l'objet
ggplot(aes(x=XXX, fill=XXX)) +
# définition du type de graphe
geom_density(aes(color=XXX), alpha=0.4) +
# réglage des couleurs de légende
scale_color_manual(values=c("XXX", "XXX")) +
scale_fill_manual(values=c("XXX", "XXX")) +
# choix de l'affichage de la légende
guides(color="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("ARNm expression") +
ylab("Density") +
labs(title="Comparison of gene expression density across different types of cancer",
fill="Cancer type")
# affichage du graphe
print(comparaison)Afficher la Solution
# création de la base de données
data <- data.frame(value=c(TCGA.PAM50_genefu_hg18$Her2,
TCGA.PAM50_genefu_hg18$Normal),
type=c(rep("Her2", 50),
rep("Normal", 50)))
# création du graphe
comparaison <- data %>%
# création de l'objet
ggplot(aes(x=value, fill=type)) +
# définition du type de graphe
geom_density(aes(color=type), alpha=0.4) +
# réglage des couleurs de légende
scale_color_manual(values=c("darkred", "darksalmon")) +
scale_fill_manual(values=c("darkred", "darksalmon")) +
# choix de l'affichage de la légende
guides(color="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("ARNm expression") +
ylab("Density") +
labs(title="Comparison of gene expression density across different types of cancer",
fill="Cancer type")
# affichage du graphe
print(comparaison)
La densité de l’expression de l’ARNmessager de ce type de cancer est beaucoup plus dispersée que celle du type de cancer Her2 et se mélange ainsi facilement dans les profils génomiques des autres types de cancer. Ce graphique met alors en évidence la difficulté à distinguer le cancer de type Normal des autres types de cancer.
Comparaison de toutes les densités
À vous de jouer
Pour mieux visualiser les différences entre tous les types de cancer, représentez la comparaison de la densité de l’ensemble des types de cancer.
Il est important de préciser au graphe que la variable doit être considéré comme une densité à l’aide de la commande aes(y=after_stat(density)).
Afficher la Solution
# création de la base de données
data <- data.frame(type=c(rep("LumA", 50),
rep("LumB", 50),
rep("Basal", 50),
rep("Her2", 50),
rep("Normal", 50)),
subtype=rep(TCGA.PAM50_genefu_hg18$chr, each=50),
value=c(TCGA.PAM50_genefu_hg18$LumA,
TCGA.PAM50_genefu_hg18$LumB,
TCGA.PAM50_genefu_hg18$Basal,
TCGA.PAM50_genefu_hg18$Her2,
TCGA.PAM50_genefu_hg18$Normal))
# densités des différents types de cancer
density <- data %>%
# création de l'objet
ggplot(aes(x=value, fill=type)) +
# définition du type de graphe
geom_density(aes(y=after_stat(density),color=type),
linewidth=0.5, alpha=0.4, fill=NA) +
# réglage des couleurs de légende
scale_color_manual(values=c("#69b3a2",
"darkred",
"darkblue",
"pink",
"darksalmon")) +
scale_fill_manual(values=c("#69b3a2",
"darkred",
"darkblue",
"pink",
"darksalmon")) +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("ARNm expression") +
ylab("Density") +
labs(title="Comparison of gene expression density across different types of cancer",
color="Cancer type")
# affichage du graphe
print(density)
Nous remarquons à nouveau une densité d’expression de l’ARNm très étalée pour le type de cancer du sein Normal. De plus, les densités pour les types de cancer Basal et LumB sont assez similaires : Basal présente un pic légèrement plus haut et et une courbe plus concentrée autour de 0.5. Globalement, chacun des types de cancer possède une expression d’ARNm concentrée à différentes valeurs pour les gènes étudiés.
POUR ALLER PLUS LOIN…
Ajoutez les histogrammes respectifs des différents types de cancer sur le graphe précédent.
Afficher la Solution
# histogrammes et densités des différents types de cancer
density <- data %>%
# créationd de l'objet
ggplot(aes(x=value, fill =type)) +
# définition du premier type de graphe
geom_histogram(aes(y=after_stat(density)),
color="#e9ecef",
alpha=0.3,
binwidth=0.5,
position="dodge") +
# définition du second type de graphe
geom_density(aes(y=after_stat(density), color=type), linewidth=0.5, fill=NA) +
# réglage des couleurs de légende
scale_color_manual(values=c("#69b3a2", "darkred", "darkblue", "pink", "darksalmon")) +
scale_fill_manual(values=c("#69b3a2", "darkred", "darkblue", "pink", "darksalmon")) +
# choix de l'affichage de la légende
guides(color="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("ARNm expression") +
ylab("Density") +
labs(title="Comparison of gene expression density across different types of cancer",
fill="Cancer type")
# affichage du graphe
print(density)
geom_histogram
Cette commande du package ggplot2 permet de créer des histogrammes sur R. De nombreuses options sont disponibles comme color, alpha (règle la transparence du remplissage), binwidth (règle la largeur des intervalles) et position = "identity", "dodge" (règle la position des bandes).
Comparaison des densités pour chacun des chromosomes
À vous de jouer
En appliquant les compétences acquises jusqu’ici, remplissez les trous afin de représenter la densité et l’histogramme d’expression de l’ARNmessager des gènes pour le type de cancer Her2 pour chacun des chromosomes.
# création de la base de données
data <- data.frame(type=XXX,
value=XXX)
# densité pour chaque chromosome
density <- data %>%
# création de l'objet
ggplot(aes(x=XXX, bins=30, color=XXX, fill=XXX)) +
# définition du premier type de graphe
geom_histogram(alpha=0.6, position="identity") +
# définition du second type de graphe
geom_density(aes(y=after_stat(ndensity)),
alpha=0.5,
fill="grey",
color="black",
linewidth=0.25) +
# séparation du graphe en sous parties
facet_wrap(~XXX, scales="XXX", nrow=3) +
# réglage des couleurs de légende
scale_fill_viridis(discrete=FALSE) +
scale_color_viridis(discrete=FALSE) +
# choix de l'affichage de la légende
guides(color="none", fill="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("Chromosome") +
ylab("Frequency") +
labs(title="Gene mRNA expression densities for Her2 type of cancer on each of the chromosomes")
# affichage du graphe
print(density)Afficher la Solution
# création de la base de données
data <- data.frame(type=TCGA.PAM50_genefu_hg18$chr,
value=TCGA.PAM50_genefu_hg18$Her2)
# densité pour chaque chromosome
density <- data %>%
# création de l'objet
ggplot(aes(x=value, bins=30, color=type, fill=type)) +
# définition du premier type de graphe
geom_histogram(alpha=0.6, position="identity") +
# définition du second type de graphe
geom_density(aes(y=after_stat(ndensity)),
alpha=0.5,
fill="grey",
color="black",
linewidth=0.25) +
# séparation du graphe en sous parties
facet_wrap(~type, scales="fixed", nrow=3) +
# réglage des couleurs de légende
scale_fill_viridis(discrete=FALSE) +
scale_color_viridis(discrete=FALSE) +
# choix de l'affichage de la légende
guides(color="none", fill="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("Chromosome") +
ylab("Frequency") +
labs(title="Gene mRNA expression densities for Her2 type of cancer on each of the chromosomes")
# affichage du graphe
print(density)
Notons que la densité ne peut pas être représentée avec un seul point ce qui explique l’absence de courbe pour les chromosomes 3, 14 et 22. Quant à lui, l’histogramme représente le nombre de gènes ayant une valeur de test à un endroit donné : par exemple, pour le chromosome 16, il y a deux gènes qui ont une valeur de test environ égale à 0.21 pour le type de cancer Her2.
Pour la plupart des chromosomes, la densité observée est bimodale. Quelques chromosomes possèdent une densité différente. En effet, le chromosome 1 possède un grand pic sur l’histogramme autour de la valeur de test 0. Le chromosome 17 est aussi atypique puisque la densité associée est assez étalée dû aux 4 pics éloignés présents de l’histogramme.
POUR ALLER PLUS LOIN…
De la même façon que le graphe précédent, remplissez les trous pour représenter graphiquement la densité et les histogrammes d’expression de l’ARNm des gènes pour chaque type de cancer sur chacun des chromosomes.
# création de la base de données pour le graphe
data <- data.frame(type=TCGA.PAM50_genefu_hg18$XXX,
subtype=c(rep("LumA", 50),
rep("LumB", 50),
rep("Basal", 50),
rep("Her2", 50),
rep("Normal", 50)),
value=c(XXX,
XXX,
XXX,
XXX,
XXX))
# densité pour chaque chromosome
density <- data %>%
# création de l'objet
ggplot(aes(x=XXX, fill=XXX, color=XXX)) +
# définition du premier type de graphe
geom_histogram(alpha=0.6, position="dodge") +
# définition du second type de graphe
geom_density(aes(y=after_stat(ndensity), color=XXX),
alpha=0.25,
fill="grey",
linewidth=0.25) +
# séparation du graphe en sous parties
facet_wrap(~XXX, scales="XXX", nrow=3) +
# réglage des couleurs de légende
scale_color_manual(values=c("#69b3a2",
"darkred",
"darkblue",
"pink",
"darksalmon")) +
scale_fill_manual(values=c("#69b3a2",
"darkred",
"darkblue",
"pink",
"darksalmon")) +
# choix de l'affichage de la légende
guides(color="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("Chromosome") +
ylab("Frequency") +
labs(title="Gene mRNA expression densities for each type of cancer on each of the chromosomes",
fill="Cancer type")
# affichage du graphe
print(density)Afficher la Solution
# création de la base de données pour le graphe
data <- data.frame(type=TCGA.PAM50_genefu_hg18$chr,
subtype=c(rep("LumA", 50),
rep("LumB", 50),
rep("Basal", 50),
rep("Her2", 50),
rep("Normal", 50)),
value=c(TCGA.PAM50_genefu_hg18$LumA,
TCGA.PAM50_genefu_hg18$LumB,
TCGA.PAM50_genefu_hg18$Basal,
TCGA.PAM50_genefu_hg18$Her2,
TCGA.PAM50_genefu_hg18$Normal))
# densité pour chaque chromosome
density <- data %>%
# création de l'objet
ggplot(aes(x=value, fill=subtype, color=subtype)) +
# définition du premier type de graphe
geom_histogram(alpha=0.6, position="dodge") +
# définition du second type de graphe
geom_density(aes(y=after_stat(ndensity), color=subtype),
alpha=0.25,
fill="grey",
linewidth=0.25) +
# séparation du graphe en sous parties
facet_wrap(~type, scales="fixed", nrow=3) +
# réglage des couleurs de légende
scale_color_manual(values=c("#69b3a2",
"darkred",
"darkblue",
"pink",
"darksalmon")) +
scale_fill_manual(values=c("#69b3a2",
"darkred",
"darkblue",
"pink",
"darksalmon")) +
# choix de l'affichage de la légende
guides(color="none") +
# personnalisation des éléments de texte
theme(strip.text.x=element_text(size=8),
plot.title=element_text(size=10, hjust=0.5, face="bold")) +
# titres
xlab("Chromosome") +
ylab("Frequency") +
labs(title="Gene mRNA expression densities for each type of cancer on each of the chromosomes",
fill="Cancer type")
# affichage du graphe
print(density)
En comparant l’expression de l’ARNm pour l’ensemble des types de cancer, nous remarquons que les conclusions faites pour le type de cancer Her2 peuvent être utilisées à nouveau. En effet, les densités associées sont généralement bimodales sauf piur quelques chromosomes. De plus, nous identifions une expression d’ARNm plus concentrée sur les chromosomes 8, 16 et 17 pour le type de cancer Her2 ainsi que sur les chromosomes 7, 11 et 12 pour le type de cancer LumA. Ce grahe permet de mettre en lumière les différentes expressions d’ARNm de mêmes gènes pour les différents types de cancer.
Création de HEATMAP
Une heatmap (ou carte thermique) est une méthode de visualisation qui permet de représenter des données de grande dimension. C’est une représentation sous forme de grille de cellules colorées. Cette méthode peut être utilisée pour visualiser les observations, les corrélations ou des valeurs manquantes par exemple.
Création d’une heatmap intéractive
L’objectif de cette partie est de réaliser des heatmap intéractives représentant les observations.
La librairie heatmaply dont nous allons nous servir est incluse dans l’écosystème de plotly (plotly) qui est un outil permettant d’inclure de l’interactivité dans les graphiques.
Dans notre exemple de heatmap, vous pourrez notament connaître la valeur d’une cellule en faisant simplement glisser votre souris dessus, zoommer sur les cellules et se déplacer dans la figure.
Pour cette première visualisation, il vous faut donc charger la librairie heatmaply et faire un choix de gradient de couleurs à l’aide des commandes suivantes :
library(heatmaply)
gradient_col <- ggplot2::scale_fill_gradient2(
low="blue", high="red", midpoint=0, limits=c(-8, 8))
print(gradient_col)<ScaleContinuous>
Range:
Limits: -8 -- 8Vous pouvez trouver ci-joint deux sites donnant des exemples d’utilisation et de l’aide sur la librairie heatmaply :
Dans un premier temps, nous allons sélectionner le chromosome que l’on souhaite étudier.
Prenons ici le chromosome numéro 1.
# on sélectionne le chromosome souhaité dans notre dataframe
chr1 <- subset(gene.exp, gene.exp$chr==1)
# on sélectionne les colonnes qui nous intéressent
mydata_1 <- chr1[, 3:18]
# nom des lignes
rownames(mydata_1) <- mydata_1[, 1]
# on supprime la première colonne qui ne contient que des 1
mydata_1 <- mydata_1[, -1]Création de la heatmap à l’aide de heatmaply :
heatmaply(mydata_1, scale_fill_gradient_fun=gradient_col,
Colv=NA, Rowv=NA, scale='none', limits=c(-8, 8),
xlab="Individus", ylab="Gènes",
main="Expression des gènes du chromosome 1")La heatmap ci-dessus nous donne visuellement une idée de l’expression de chaque gène du chromosome 1. On observe que dans la majorité des cas le gène CLCA2 est sur-exprimé (les cellules sont rouges), et les gènes PDZK1 DLANI1 sont sous-exprimés (les cellules sont bleutées).
Vous pouvez choisir d’afficher les corrélations avec la fonction heatmaply_cor() ou les valeurs manquantes avec la fonction heatmaply_na().
A vous de jouer
En vous basant sur l’exemple ci-dessus, créez une heatmap représentant l’expression des gènes du chromosome 17.
Voici une idée du résultat que vous devez obtenir :
Afficher la Solution
chr17 <- subset(gene.exp, gene.exp$chr == 17 )
mydata_17 <- chr17[,3:18]
rownames(mydata_17) <- mydata_17[,1]
mydata_17 <- mydata_17[,-1]
heatmaply(mydata_17, scale_fill_gradient_fun=gradient_col,
Colv=NA, Rowv=NA, scale='none', limits=c(-7, 8),
xlab="Individus",ylab="Gènes",
main="Expression des gènes du chromosome 17")Création d’une heatmap circulaire
L’objectif de cette seconde partie sur les heatmap est de réaliser des heatmap circulaires avec dendrogrammes.
Note
La méthode de clustering appliquée est hclust soit une méthode de classiification hierarchique, la distance utilisée par défaut est celle de Ward. Vous pouvez trouver des informations sur hclustau Lien suivant et sur la classification hierarchique ici et ici.
Commencez par charger les packages nécessaires :
library(circlize)
library(ComplexHeatmap)On effectue une petite modification de notre dataframe en supprimant une colonne qui ne nous sera pas utile par la suite.
# suppression de la colonne position
gene.exp2 <- gene.exp[, -2]Pour pouvoir réaliser notre heatmap, il nous faut tout d’abord convertir notre dataframe en matrice :
mat <- as.matrix(gene.exp2, labels=TRUE) # on transforme le dataframe en matrice
# On définit les noms des colonnes et des lignes
row_names <- as.vector(gene.exp2$NAME)
col_names <- c("chr", "NAME", "TCGA.A2.A04W.01A", "TCGA.A2.A0D1.01A",
"TCGA.A2.A0EQ.01A", "TCGA.A8.A081.01A", "TCGA.A8.A08J.01A",
"TCGA.A8.A09X.01A", "TCGA.AN.A0FV.01A", "TCGA.AO.A0J2.01A",
"TCGA.AO.A12D.01A", "TCGA.B6.A0RS.01A", "TCGA.BH.A0EE.01A",
"TCGA.C8.A12L.01A", "TCGA.C8.A12P.01A", "TCGA.D8.A13Z.01A",
"TCGA.E2.A1B0.01A")
rownames_factor <- factor(row_names)
# On met notre matrice sous format "numeric"
mat_numeric <- matrix(as.numeric(mat), ncol = ncol(mat))
rownames(mat_numeric) <- row_names
colnames(mat_numeric) <- col_names
# Vecteur des chromosomes :
split <- as.vector(gene.exp2$chr)
split <- factor(split)Le vecteur splitdoit contenir les différents groupes/catégories si vous souhaitez divisez votre heatmap par groupes/catégories.
Création des heatmap circulaires :
- Sur cette première heatmap, nous avons affiché tous les chromosomes ainsi que tous les gènes de manière simple.
## Affichage des heatmaps
circos.clear()
circos.par(start.degree=90, gap.degree=3)
col_fun1 = colorRamp2(c(-6, 0, 6), c("blue", "white", "red"))
# heatmap simple
circos.heatmap(mat_numeric[, !(colnames(mat_numeric) %in% c("chr", "NAME"))],
split=split,
col=col_fun1)
- Sur cette seconde heatmap, nous avons encadré chaque chromosome et affiché leur numéro.
# heatmap avec les chromosomes associés
circos.clear()
circos.heatmap(mat_numeric[, !(colnames(mat_numeric) %in% c("chr", "NAME"))],
split=split, col=col_fun1, track.height=0.4,
bg.border="green", bg.lwd=2, bg.lty=2,
show.sector.labels=TRUE)
- Sur cette troisième heatmap, nous avons rajouté le nom de chaque gène.
# heatmap avec les chromosomes et les noms des gènes
circos.clear()
circos.heatmap(mat_numeric[, !(colnames(mat_numeric) %in% c("chr", "NAME"))],
split=split, col=col_fun1, track.height=0.4,
bg.border="green", bg.lwd=2, bg.lty=2,
show.sector.labels=TRUE,
rownames.side="outside", rownames.cex=0.4)
- Sur cette quatrième heatmap, nous avons rajouté les dendrogrammes ainsi qu’une légende de couleur.
# heatmap avec les dendrogrammes
circos.clear()
circos.heatmap(mat_numeric[, !(colnames(mat_numeric) %in% c("chr", "NAME"))],
split=split, col=col_fun1, track.height=0.4,
bg.border="green", bg.lwd=2, bg.lty=2,
show.sector.labels=TRUE,
dend.side="outside")
lgd = Legend(title="values", col_fun=col_fun1)
grid.draw(lgd)
Le nombre de chromosomes étant important, il est compliqué de voir en détail les gènes. L’exercice suivant consiste donc à sélectionner certains chromosomes et à afficher la heatmap circulaire correspondante.
Il est important de toujours appeler la commande circos.clear() avant chaque nouvelle heatmap.
A vous de jouer
Sélectionnez les chromosomes 1, 3, 6, 11, 17 et 19.
x <- c("XXX") # vecteur des chromosomes que l'on souhaite sélectionner
data_chr <- gene.exp2[gene.exp2$chr %in% x, ] # dataframe avec nos chromosomes sélectionnés
# Convertion en matrice
mat_chr <- as.matrix(data_chr, labels=TRUE)
mat_chr_numeric <- matrix(as.numeric(mat_chr), ncol = ncol(mat_chr))
# On définit les noms des colonnes et des lignes
row_names_chr <- as.vector(data_chr$NAME)
col_names_chr <- c("chr", "NAME", "TCGA.A2.A04W.01A", "TCGA.A2.A0D1.01A",
"TCGA.A2.A0EQ.01A", "TCGA.A8.A081.01A", "TCGA.A8.A08J.01A",
"TCGA.A8.A09X.01A", "TCGA.AN.A0FV.01A", "TCGA.AO.A0J2.01A",
"TCGA.AO.A12D.01A", "TCGA.B6.A0RS.01A", "TCGA.BH.A0EE.01A",
"TCGA.C8.A12L.01A", "TCGA.C8.A12P.01A", "TCGA.D8.A13Z.01A",
"TCGA.E2.A1B0.01A")
rownames(mat_chr_numeric) <- row_names_chr
colnames(mat_chr_numeric) <- col_names_chr
# Vecteur des chromosomes :
split2 <- as.vector(data_chr$chr)
split2 <- factor(split2)Réaliser maintenant une heatmap affichant les noms des gènes ainsi que les dendrogrammes.
Voici une idée du résultat que vous devez obtenir :
Afficher la Solution
x <- c(1, 3, 6, 11, 17, 19)
data_chr <- gene.exp2[gene.exp2$chr %in% x, ]
mat_chr <- as.matrix(data_chr, labels=TRUE)
mat_chr_numeric <- matrix(as.numeric(mat_chr), ncol = ncol(mat_chr))
row_names_chr <- as.vector(data_chr$NAME)
col_names_chr <- c("chr", "NAME", "TCGA.A2.A04W.01A", "TCGA.A2.A0D1.01A",
"TCGA.A2.A0EQ.01A", "TCGA.A8.A081.01A", "TCGA.A8.A08J.01A",
"TCGA.A8.A09X.01A", "TCGA.AN.A0FV.01A", "TCGA.AO.A0J2.01A",
"TCGA.AO.A12D.01A", "TCGA.B6.A0RS.01A", "TCGA.BH.A0EE.01A",
"TCGA.C8.A12L.01A", "TCGA.C8.A12P.01A", "TCGA.D8.A13Z.01A",
"TCGA.E2.A1B0.01A")
rownames(mat_chr_numeric) <- row_names_chr
colnames(mat_chr_numeric) <- col_names_chr
split2 <- as.vector(data_chr$chr)
split2 <- factor(split2)
circos.clear()
circos.heatmap(mat_chr_numeric[, !(colnames(mat_chr_numeric) %in% c("chr", "NAME"))],
split=split2, col=col_fun1, track.height=0.4,
bg.border="green", bg.lwd = 2, bg.lty=2,
show.sector.labels=TRUE, dend.side="outside",
rownames.side="inside", rownames.cex=0.35)
lgd = Legend(title="values", col_fun=col_fun1)
grid.draw(lgd)
Vous pouvez consulter le lien suivant si vous souhaitez en apprendre plus sur les heatmap circulaires : heatmap circulaire
Visualisation avec Circos
A présent, nous allons réaliser une visualisation circulaire à l’aide du package OmicCircos et de la fonction circos().
Vous pouvez trouver de l’aide au package OmicCircos aux liens suivants :
colors <- rainbow(10, alpha = 0.5) # choix des couleurs
# Si vous souhaitez enregistrer votre image sous forme de pdf
#pdf("visucomp_histo.pdf", 8,8) # nom de l'image
#par(mar=c(2, 2, 2, 2)) # definition de la fenetre graphique
plot(c(1,800), c(1,800), type="n", axes=FALSE, xlab="", ylab="", main="")
zoom <- c(1, 22, 939245.5, 154143883, 0, 180)
###### Moitié droite du cercle ##########
# mettre les chromosomes et l'emplacement des gènes
circos(type="chr", R=400, cir="hg18", W=4,
print.chr.lab=TRUE, scale=TRUE, zoom=zoom)
# création de la heatmap
circos(type="heatmap2", R=300, cir="hg18", W=100,
mapping=gene.exp, col.v=4, cluster=FALSE,
col.bar=TRUE, lwd=0.01, zoom=zoom)
# nombre de copies (gain_rouge, perte_bleu)
circos(type="ml3", R=220, cir="hg18", W=80, mapping=cnv,
col.v=4, B=FALSE, lwd=1, cutoff=0,
zoom=zoom)
# création de l'histogramme
circos(type="h", R=140, cir="hg18", W=80, mapping=cnv,
col.v=4, B=TRUE, lwd=1, col=colors[1],
zoom=zoom)
# Mise en valeur des chromosomes 11 et 17
the.col1 = rainbow(10, alpha=0.5)[1]
highlight = c(140, 400, 11, 282412.5, 11, 133770314.5, the.col1, the.col1)
circos(type="hl", R=110, cir="hg18", W=40,
mapping=highlight, lwd=2, zoom=zoom)
the.col2 = rainbow(10, alpha=0.5)[6];
highlight = c(140, 400, 17, 739525, 17, 78385909, the.col2, the.col2)
circos(type="hl", R=110, cir="hg18", W=40,
mapping=highlight, lwd=2, zoom=zoom)
highlight.link1 = c(400, 400, 140, 376.8544, 384.0021,
450, 540.5)
circos(type="highlight.link", cir="hg18",
mapping=highlight.link1, col=the.col1, lwd=1)
highlight.link2 = c(400, 400, 140, 419.1154, 423.3032,
543, 627)
circos(type="highlight.link", cir="hg18",
mapping=highlight.link2, col=the.col2, lwd=1)
####### Moitié gauche du cercle #######
# Chromosome 11
zoom = c(11, 11, 282412.5, 133770314.5, 180, 270)
circos(type = "chr",R = 400, cir = "hg18", W = 4, print.chr.lab = TRUE,
scale = TRUE, zoom = zoom)
circos(type = "heatmap2", R = 300, cir = "hg18", W = 100,
mapping = gene.exp, col.v = 4, cluster = FALSE, lwd = 0.01,
zoom = zoom)
circos(type = "ml3", R = 220, cir = "hg18", W = 80, mapping = cnv,
col.v = 4, B = FALSE, lwd = 1, cutoff = 0,
zoom = zoom)
circos(type = "h", R = 140, cir = "hg18", W = 80, mapping = cnv,
col.v = 4, B = TRUE, lwd = 1, col = colors[1],
zoom = zoom)
# Chromosome 17
zoom = c(17, 17, 739525, 78385909, 274, 356)
circos(type = "chr", R = 400, cir = "hg18", W = 4, print.chr.lab = TRUE,
scale = TRUE, zoom = zoom)
circos(type = "heatmap2", R = 300, cir = "hg18", W = 100,
mapping = gene.exp, col.v = 4, cluster = FALSE, lwd = 0.01,
zoom = zoom)
circos(type = "ml3", R = 220, cir = "hg18", W = 80, mapping = cnv,
col.v = 4, B = FALSE, lwd = 1, cutoff = 0,
zoom = zoom)
circos(type = "h", R = 140, cir = "hg18", W = 80, mapping = cnv,
col.v = 4, B = TRUE, lwd = 1, col = colors[1],
zoom = zoom)
De l’extérieur vers l’intérieur, vous pouvez trouver :
le numéro du chromosome et la position des gènes
une heatmap de l’expression des gènes
le nombre de copies de chaque gènes par rapport à la normale (rouge : gain, bleu : perte)
des histogrammes pour observer la dispersion de la loi de distribution
Lexique
hg18 : La représentation du génome humain a connu plusieurs versions avant d’aboutir à la représentation la plus optimale. La version 18 correspond à l’avant dernière version.
Her2 : Il s’agit d’une protéine naturellement présente dans l’organisme : un récepteur transmembranaire impliqué dans la régulation de la prolifération cellulaire.
LumA : Le cancer du sein Luminal A est l’un des sous-types luminaux et est généralement associé à des caractéristiques moins agressives que Luminal B.
LumB : Ce type de cancer du sein est associé à un grade plus élevé, des taux de prolifération accrus, et un pronostic global plus défavorable.
Basal : Il s’agit d’un sous-type du cancer du sein négatif pour les récepteurs d’oestrogène (ER), les récepteurs de progestérone (PR) et le récepteur 2 du facteur de croissance épidermique humain (HER2). On le désigne souvent comme un cancer du sein triple négatif (TNBC).
FDR (First Division Restitution) : Ceci réfère à l’évènement où l’une des cellules filles produites lors de la première division méiotique conserve les deux chromatides d’un chromosome, sans subir la séparation normale en deux cellules distinctes. Cela aboutit à un gamète avec un ensemble complet de chromosomes, plutôt que la moitié attendue.
Bonferroni : Cette correction statistique permet d’ajuster le niveau de confiance de chaque intervalle pour que le niveau de confiance simultané ne soit pas erroné.
Graphiques supplémentaires
Voir le diagramme en barre
data(TCGA.PAM50_genefu_hg18)
df<-TCGA.PAM50_genefu_hg18
# Pour voir mieux, j'ai affiché que les 10 premières lignes du dataset. Si vous voulez voir plus, vous pouvez mettre en commentaire la ligne suivante.
df<-df[1:10,]
df_long <- df %>%
gather(key = "Subtype", value = "Expression", -chr, -po, -Gene)
ggplot(df_long, aes(x = Gene, y = Expression, fill = Subtype)) +
geom_bar(stat = "identity") +
labs(title = "Expression génique des gènes dans les sous-catégories de cancer du sein",
x = "Gène", y = "Expression génique") +
theme_minimal()